inline	
grass_kernels::grass_kernels(const field<mat> in_grass_points)
:grass_points(in_grass_points), N_points (in_grass_points.n_rows)
{
  
}

inline 
mat 
grass_kernels::projection()	// Projection Kernel - Sareh's paper (Clustering on Grassmann....)
{
  
  cout << "Calculating Projection Kernel using Grassmann Points: " << endl;  
  mat K = zeros(N_points, N_points);
  
  mat tmp;
  for (uword i = 0; i < N_points ;  ++i)
    for (uword j = 0; j < N_points ;  ++j)
    {
      /*
       * Dont do it in this way: Insane, slower :(
       * tmp = grass_points(i,0);
      mat tmp_i = tmp*tmp.t(); 
      tmp = grass_points(j,0);
      mat tmp_j = tmp*tmp.t(); 
      K(i,j) = trace(tmp_i*tmp_j);
      */
     tmp = grass_points(i,0).t()*grass_points(j,0);
     K(i,j) = pow( norm(tmp,"fro"), 2);
     
      
    }
    
    
    //K.print("K no fast");
    //K.save("Knofast.dat", raw_ascii);
    //projection2();
    return K;
}

inline 
mat 
grass_kernels::projection2()	// Projection Kernel - Sareh's paper (Clustering on Grassmann....)
{
  
  cout << "Calculating Projection Kernel 2-- Fast!!!!!!!!: " << endl;  
  mat K = zeros(N_points, N_points);
  mat K2 = zeros(N_points, N_points);// for comparison
  
  mat tmp;
  cout << "i = " << endl;
  for (uword i = 0; i < N_points ;  ++i)
  {
    cout << i << "... ";
    for (uword j = i; j < N_points ;  ++j)
    {
      /*
       * Dont do it in this way: Insane, slower :(
      tmp = grass_points(i,0);
      mat tmp_i = tmp*tmp.t(); 
      tmp = grass_points(j,0);
      mat tmp_j = tmp*tmp.t(); 
      K(i,j) = trace(tmp_i*tmp_j);
      */
     tmp = grass_points(i,0).t()*grass_points(j,0);
     K(i,j) = pow( norm(tmp,"fro"), 2 );
      
    }
  }
    
    //K.print("K before");
    //K.save("Kafter.dat", raw_ascii);
    K = symmatu(K);
    //K.print("K after");
    //K.save("Kfast.dat", raw_ascii);
    return K;
}



///Binet-Cauchy Kernel
inline 
mat 
grass_kernels::binet_cauchy()	// Binet-Cauchy Kernel - Mehrtash's paper (Sparse Representation over Grassmann Manifolds...)
{
  cout << "Calculating Binet-Cauchy Kernel using Grassmann Points: " << endl;  
  mat K = zeros(N_points, N_points);
  mat tmp;
  
  for (uword i = 0; i < N_points ;  ++i)
    for (uword j = 0; j < N_points ;  ++j)
    {
      mat mul1 = grass_points(i,0).t()*grass_points(j,0);
      mat mul2 = grass_points(j,0).t()*grass_points(i,0);
      
      //mul1.save("mul1.mat",raw_ascii);
      //mul2.save("mul2.mat",raw_ascii);
      //grass_points(i,0).save("grass_points_i.mat",raw_ascii);
      //grass_points(j,0).save("grass_points_j.mat",raw_ascii);
      double val;
      double sign;
      log_det(val, sign, mul1*mul2);
      K(i,j) = exp(val)*sign;
      //cout << "K(i,j) A: " << K(i,j)  << endl;
      //getchar();
      //cout << "K(i,j) B: " << det(mul1*mul2)  << endl;
      //getchar();
      //cout << i << " - " << j << endl;
    }
    //K.save("Kbc.mat",raw_ascii);
    //K.print("Kbc");
    return K;
}


///Canonical Correlation Pseudo-Kernel
inline 
mat 
grass_kernels::can_corr()	// Canonical Correlation Pseudo-Kernel - - Mehrtash's paper (Sparse Representation over Grassmann Manifolds...)
{
  cout << "Calculating Canonical Correlation Pseudo-Kernel using Grassmann Points: " << endl;  
  mat K = zeros(N_points, N_points);
  
  for (uword i = 0; i < N_points ;  ++i)
    for (uword j = 0; j < N_points ;  ++j)
    {
      mat Y;
      Y = grass_points(i,0).t()*grass_points(j,0); // as per (Eq. 4.2) in Jihun Hamm's Thesis (2008)
      mat U;
      vec s;
      mat V;
      svd(U,s,V,Y);  // s is organised in a decreasing order. 
      
      //Kcc(X,Y) = cos(theta1). Theta1 is the smallest principal angle. 
      //We need the first value of s, as s = ( cos(theta1) ... cos(thetam) ))
      
      
      //s.print("s:");
      
      K(i,j) = s(0);// cos(theta)
      //cout << "K(i,j) " << K(i,j)  << endl;
      //cout << i << " - " << j << endl;
    }
    //K.print("Canonical Correlation Pseudo-Kernel");
    return K;
    
    
    
}